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Resumen 


El drenaje subterráneo es utilizado para eliminar excedentes de agua en la zona 
radical y suelos salinos para lixiviar las sales. La dinámica del agua es estudiada 
con la ecuación de Boussinesq; sus soluciones analíticas son obtenidas asumiendo 
que la transmisibilidad del acuífero y la porosidad drenable son constantes, y que 
la superficie libre se abate de manera instantánea sobre los drenes. La solución en 
el caso general requiere de soluciones numéricas. En la literatura se ha demostrado 
que la condición de frontera en los drenes es una condición de radiación fractal y 
que la porosidad drenable es variable y relacionada con la curva de retención de 
humedad, y ha sido resuelta con el método del elemento finito, que en un esquema 
unidimensional puede hacerse equivalente al método de diferencias finitas. 
Aquí se propone una solución en diferencias finitas de la ecuación diferencial, 
considerando la porosidad drenable variable y la condición de radiación fractal. 
El esquema en diferencias finitas propuesto ha resultado en dos formulaciones: 
en una aparecen de manera explícita la carga y porosidad drenable, variables 
ligadas con una relación funcional, que se ha denominado esquema mixto; en 
la otra aparece sólo la carga hidráulica, denominada esquema en carga. Los dos 
esquemas coinciden cuando la porosidad drenable es independiente de la carga. 
Los esquemas han sido validados con una solución analítica lineal, y para la 
no linealidad se ha mostrado que es estable y concisa. La solución numérica es 
útil para la caracterización hidrodinámica del suelo a través de una modelación 
inversa, y para un mejor diseño de los sistemas de drenaje agrícola subterráneo, 
ya que las hipótesis consideradas en las soluciones clásicas han sido eliminadas. 


Palabras clave: formulación mixta, formulación en carga, curva de retención, 
modelación inversa. 


Introducción 


Los sistemas de drenaje subterráneos son 
ampliamente utilizados en la agricultura para 
eliminar excedentes de agua en la zona radical 
de las plantas y para lixiviar las sales del perfil 
de los suelos. La dinámica del agua en estos 
sistemas ha sido estudiada aceptando la validez 
de la ley de Darcy (1956) y en función de la escala 
de estudio se pueden utilizar dos ecuaciones 
diferenciales. La ecuación de Richards (1931) 


—que resulta de la aplicación del principio de 
conservación de la masa en el flujo del agua en 
un volumen elemental de medio poroso y de la 
ley de Darcy— permite considerar la geometría 
de los drenes en las condiciones de frontera; 
sin embargo, la simulación de la dinámica 
del agua con las soluciones numéricas bi o 
tridimensionales puede ser ardua (Zavala et al., 
2003). La ecuación de Boussinesq (1904) de los 
acuíferos libres —que resulta de la aplicación 
del principio de conservación de la masa en una 
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columna elemental de medio poroso y de la 
propia ley de Darcy— pondera las propiedades 
del suelo y el sistema en la vertical, y es a lo 
más una ecuación bidimensional; el acuífero se 
modela en planta y la geometría de los drenes 
es introducida como líneas matemáticas oO 
como puntos en un análisis bidimensional o 
unidimensional, respectivamente. 

La ecuación de Boussinesq unidimensional 
ha sido una de las bases para construir 
soluciones analíticas aproximadas de la 
dinámica del agua en un sistema de drenaje 
tanto en régimen permanente como transitorio 
(e.g. Hooghoudt, 1940; Dumm, 1954), y que 
son utilizadas en el diseño de los sistemas. En 
la deducción de la ecuación de Glover-Dumm 
para el régimen transitorio se asume que la 
transmisibilidad del acuífero y la porosidad 
drenable son constantes, y que la superficie 
libre se abate de manera instantánea sobre 
los drenes. Dado que los tres supuestos no 
representan adecuadamente las condiciones 
reales, la solución en cuestión puede ser de 
aplicabilidad limitada. Sin embargo, consi- 
derando condiciones más representativas, 
conduce a dificultades analíticas, razón por 
la cual es necesaria la utilización de métodos 
numéricos para construir soluciones de la 
ecuación de Boussinesq. 

En una línea de investigación, Zavala et al. 
(2004, 2007) analizan detalladamente el tipo 
de condición de frontera representativa de 
las condiciones reales; los autores, basados 
en los conceptos de la geometría fractal y en 
experiencias de drenaje, recomiendan una 
condición de radiación fractal, la cual incluye 
la radiación lineal utilizada por Fuentes et al. 
(1997). En cuanto a la porosidad drenable, 
Fuentes et al. (2009), basados en los conceptos 
de lámina drenable y lámina drenada, y 
también en experiencias de drenaje, propo- 
nen una expresión analítica en la cual 
interviene la curva de retención de humedad 
de los suelos. Los autores citados utilizan el 
método del elemento finito para resolver la 
ecuación unidimensional de Boussinesq, con 
buenos resultados en cuanto a la estabilidad, 


convergencia y precisión de la solución. En 
un esquema unidimensional, el método del 
elemento finito puede hacerse equivalente 
al método de diferencias finitas (Russell y 
Wheeler, 1983). 

Existe un esquema en diferencias finitas, 
basado en el esquema de Laasonen, propuesto 
por Zataráin et al. (1998), para resolver 
Richards 
aplicada al fenómeno de la infiltración del 


numéricamente la ecuación de 
agua en los suelos, con excelentes resultados. 
Aparte de su alta precisión, estabilidad y 
convergencia, el esquema tiene la ventaja 
adicional de su naturaleza intuitiva, ya que 
está basado en un balance local de masa. Este 
esquema puede ser utilizado para resolver la 
ecuación unidimensional de Boussinesq del 
drenaje agrícola. 

El presente trabajo tiene como objetivo 
la solución numérica de la ecuación unidi- 
mensional de Boussinesq con el método de 
diferencias finitas basado en un balance local de 
masa. Se considera que la porosidad drenable 
es variable y las condiciones de frontera en los 
drenes son de radiación fractal. 


Ecuaciones de base 
Ecuación de Boussinesq 


En el estudio de la dinámica del agua en 
sistemas de drenaje con la ecuación de 
Boussinesq se asume generalmente que las 
variaciones de la carga hidráulica a lo largo 
de los tubos de drenaje (dirección y) son 
despreciables respecto a las variaciones de 
la carga en un corte transversal (dirección x). 
En esta situación, la ecuación a resolver en el 
dominio mostrado en la figura 1 es la ecuación 
de Boussinesq unidimensional, que resulta de 
la ecuación de continuidad: 


200) 2 (1) 
ot 0x 
y de la ley de Darcy: 
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0H 
=-K,=— 2 
q > (2) 
a saber: 
0H Y 0H 
H T(H R 3 
O 


donde H = HÍ(x, t) es la carga hidráulica 
contada a partir de un estrato impermeable, y 
es una función de la coordenada horizontal (x) 
y del tiempo (t); q es el flujo de Darcy o caudal 
por unidad de área; K, es la conductividad 
hidráulica a saturación; v=v(H) es la porosidad 
drenable como una función de la carga; R,, es 
el volumen de recarga en la unidad de tiempo 
por unidad de área del acuífero. 

El caudal unitario de agua Q, = Hg es pro- 
porcionado por: 


0 =Tmá ; T(B=K.H (4) 


donde T(H) es la transmisibilidad. 
La capacidad de almacenamiento está defi- 
nida por: 


0 (5) 


Figura 1. Esquema de un sistema de drenaje agrícola 


subterráneo. 


donde W = (vH) es la lámina de agua drenable. 
La igualdad u = v se da cuando la porosidad 
drenable es independiente de la carga. 


La porosidad drenable 


Siguiendo el procedimiento de Fuentes et al. 
(2009), se establece una expresión de la capa- 
cidad de almacenamiento. Si la superficie 
libre estuvo inicialmente en la posición z = H, 
la lámina drenada, cuando esta superficie 
se encuentra en la posición z = H< H, está 


definida por: 


L(H)= f(o.- Op (2) ]4z (6) 


donde 0, es el contenido de humedad a 
saturación; 0,,, (2) es el contenido de humedad 
en función de la posición o perfil de humedad. 

De acuerdo con Fragoza et al. (2003), la 
capacidad de almacenamiento de un acuífero 
libre se define como: 


1 (H)=0,-0(H-H,) (7) 


donde 0, es el contenido de humedad a 
saturación; O (H — H,) representa la evolución 
del contenido de humedad en la posición z = 
H, mientras la superficie libre desciende; z es la 
elevación de la superficie del terreno. 

La porosidad drenable se deduce a partir 
de la igualdad de las ecuaciones (5) y (7), a 
saber: 


0009 [u(jai=2 [oo (in )jañ 6 


donde H es la variable de integración. 

Para calcular la capacidad de almace- 
namiento y la porosidad drenable es necesario 
proporcionar la curva de retención de hume- 
dad del suelo. En la literatura es bastante 
común representarla con la ecuación de van 
Genuchten (1980): 
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9(w)=0,+(0, 0) 8] (9) 


donde 6, es el contenido residual de humedad; 
y, es un parámetro de escala de la presión; m y 
n son dos parámetros de forma positivos. 

La introducción de la ecuación (9) en las 
ecuaciones (7) y (8) proporciona la capacidad 
de almacenamiento y la porosidad drenable 
siguientes: 


u(H)=(0,-0,),1 Esa] 


di íÁ (+ y)" “a Jn 


Hu -Hy/lva| 


v(B)=(0, cop lval 


La porosidad drenable no tiene una forma 
analítica cerrada, pudiendo ser calculada me- 
diante integración numérica. 

Una forma cerrada puede ser construida a 
partir de la difusividad de Fujita (1952) y de la 
relación conductividad hidráulica-difusividad 
hidráulica de Parlange et al. (1982) (Fuentes et 
al., 1992), definidas por: 


KA 1-0 
D(0) = == 
a (1-00) 


e|1-$+(P-aje] 
1-00 


(12) 


K(0)=K., (13) 


donde O = (0-0)/ (0-0) es un grado efectivo 
de saturación; O. y PB son parámetros de forma 
adimensionales, tales que0d<a<1ly0<B<l; 
A, es la escala de Bouwer (1964). 

De la definición de la difusividad hidráulica 
D(0) = K(0)dy/d0, considerando la condición 0 
= 0 cuando y = 0, se deduce: 


DN Bb As 50) 


y BS Q al 
“B(=B) (1- ajo 


(14) 


donde y.=-1.. 
La porosidad drenable se obtiene de las 
ecuaciones (8) y (14): 
e(H-H;) 
(15) 


o(H)=(0.- 2 tn e 


1-00 


Se debe notar que la función 0(y) es implí- 
cita en la ecuación (14), y en consecuencia 
la función v(H). Estas funciones pueden ser 
explicitadas en función de la presión si se 
acepta QU. = P; en este caso, la conductividad 
en función de la presión corresponde a la 
ecuación de Gardner (1958) K(w) = K_ exp(wy/A.), 
ampliamente utilizada en estudios teóricos. La 
curva 0(y) correspondiente es la siguiente: 


6.-6 
6 =0. AA, 16 
Wa 


De la ecuación (7) se obtiene la capacidad de 
almacenamiento: 


(4) =(0, 0) a 022) Y 


y de la ecuación (15), la porosidad drenable: 


v(H)=(0,-0,) 
l o qe o. +0.exp[ (H-H,)/», al (18) 


1-0, +0 exp[-H, /A, ] 


La lámina drenada se obtiene de la ecuación 


(6): 
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e(5)=(6.- 09100 


lla]! 


El contenido de humedad a saturación pue- 
de ser asimilado a la porosidad total (p), la cual 
es estimada a partir de la densidad total del 
suelo seco (p,) y de la densidad de las partículas 
(p,), con la fórmula $ =1 — p,/ p,; el contenido 
de humedad residual puede ser asumido igual 


(19) 


a Cero. 
Condiciones inicial y de frontera 


La carga hidráulica contada a partir del estrato 
impermeable H(x, t) está relacionada con la 
carga h(x, t) contada a partir de los drenes, de 
acuerdo con la figura 1, por: 


H(x,t)=D,+h(x,t) (0) 
donde D, es la altura de los drenes a partir del 
estrato impermeable. 

La variación transversal de h al inicio del 
proceso de drenaje es considerada como la 
condición inicial: 

h(x,0) =h,(x) (21) 

En cuanto a las condiciones de frontera o 
condiciones en los drenes ubicados en x = 0 
y x = L, se han asumido formas diversas. La 
solución de Glover-Dumm es establecida 
asumiendo que la carga sobre el dren se abate 
totalmente de manera instantánea (Dumm, 
1954); esta condición es de tipo Dirichlet o 
de primer orden. La solución de Fuentes et 
al. (1997) (ver también Fragoza et al., 2003) se 
obtiene bajo el argumento soportado expe- 
rimentalmente de que el flujo de Darcy en 
los drenes es proporcional a la carga (q «< h); 
esta condición es de tipo radiación lineal o de 
tercer orden. Ambas condiciones de frontera 
son casos particulares de: 


(22) 


donde el signo positivo corresponde al dren 
posicionado en x = 0; mientras que el negativo, 
al posicionado en x = L. 

El flujo de radiación lineal se puede 
expresar como q = KKh/L, donde k es un 
coeficiente adimensional de conductancia de la 
interfaz suelo-dren. De la solución de Fuentes 
et al. (1997), se deduce la solución de Glover- 
Dumm cuando este coeficiente es infinito. En 
esta misma línea de investigación, Zavala et 
al. (2004, 2007) proponen una condición de 
radiación fractal: 


2s 
ko) =0. ; x=0,L (23) 


donde q, es el flujo correspondiente y depende 
de las características de la interfaz suelo-dren; 
s está definido por s = D/E, donde D es la 
dimensión fractal efectiva de la interfaz suelo- 
dren y E = 3 es la dimensión de Euclides del 
espacio físico. La relación entre s y la porosi- 
dad efectiva de la interfaz es proporcionada 
por la ecuación presentada por Fuentes et al. 
(2001): 
(1-0)+9%=1 (24) 
En un sistema de drenes paralelos a igual 
separación, el gasto de agua que fluye a través 
de la frontera por unidad de longitud de dren 
es (Fuentes et al., 1997): 


Q,(4)=2[D,+»(0,£)]4, [1(0,2/n,] (25) 


La evolución temporal de la lámina drenada 
se calcula con la siguiente expresión: 


10 Jo.(Dar (26) 


donde tes una variable de integración. 
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Solución numérica 
Esquemas numéricos 


La ecuación unidimensional de Boussinesq se 
resuelve con el método de diferencias finitas, 
adaptando el esquema numérico propuesto 
por Zataráin et al. (1998) para un problema 
similar en la escala de la ecuación de Richards. 
La adaptación a la escala de la ecuación de 
Boussinesq requiere de la discretización del 
dominio, como se muestra en la figura 2. 

Para plantear la resolución numérica de la 
ecuación (3), se introducen los parámetros de 
interpolación definidos por: 


a PEZ qa OY (27) 
Xi Eat, 


tales que0<y<1y0<0<1;1=1,2,...y¡=172,... 
son los índices para el espacio y el tiempo, 
respectivamente. 

La variable dependiente (H) en un nodo 
intermedio 1 + y para todo j se estima como: 


us [0 —y)H!+ YH, (28) 


mientras que en el tiempo intermedio / + () para 
todo ise estima como: 


H*"= (1-0)H/+ 0H" (29) 


La ecuación (1) se aplica en el tiempo En 


la derivada temporal puede ser discretizada de 
acuerdo con las dos formulaciones siguientes: 


+0 ME 5248 
A(AJ"_ ni ol] 
dt | At, 


1 ] 


Ati= ti És 


, 60) 


1 


1 
ot ¡ At; 


H ¡+0 JH ri 
9(vA) Epa (31) 


Para identificar en lo sucesivo los dos 
esquemas numéricos resultantes, el primero 
es denominado esquema mixto y el segundo 
esquema en carga, ya que en el primero 
aparecen de manera explícita la carga y el 
volumen de agua que se drena, mientras que 


Figura 2. El dominio de resolución de la ecuación de Boussinesq; (w es el factor de interpolación en el tiempo 


y y es el factor de interpolación en el espacio. 
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en el segundo aparece explícitamente sólo la 
carga; las dos formulaciones coinciden cuando 
la porosidad drenable es independiente de la 
carga y la formulación en carga no requiere de 
la integración numérica en la ecuación (11) para 
calcular la porosidad drenable. 

La discretización de la derivada espacial 
alrededor del nodo ¡-ésimo es la siguiente: 


(ag) May, 
ox | Ax; ó 


Í 1 


ax,= (1-1); 


(32) 


X1)+ 10,4 x;) 


De la ecuación (4) se obtiene el caudal 
unitario definido en los nudos intermedios: 


o + HE 1 Hp + 
a, iy TE) (83) 


1+Y 
HI 


(Hg) 


i+Y 


j+0 
Hoy 


TER =T(HEGp) 


. J40_ piro 
O H; Hi e 
FIA) y _ ñ 

XX Gx 


(Hg) 


(34) 


Las cargas en los diferentes nodos y en el 
tiempo intermedio se obtienen de la ecuación 
(09), 
ecuaciones (33) y (34), y éstas, a su vez, en la 


las cuales son introducidas en las 
ecuación (32). Luego, las ecuaciones (30) y (32) 
se llevan a la ecuación (1), y se asocian términos 
semejantes, resultando el sistema de ecuaciones 
algebraicas siguiente: 


AH BHPS+D HA =E,; ¿i=2,3,...n-1 (35) 

donde: 
OT 
A¡=- A (36) 
Ax; (x;-x;1) 
j+ pei j 
B, o| Th a) |, 0% (37) 
AX¡| XX, X=X,1 |] At 


OTÍ+o 
D.= + (38) 
Ax; (x HT x;) 


E; =R j+o y ls 0) TER Ha TE a 


11 
Ax; 


| ea Xi 
(39) 
j j+0 j+o 
o_o) TES, Tas 137 
pa 
Abi AX | A 2% 


Para el esquema en carga, ecuación (31), 
los coeficientes B, y E, deben ser redefinidos 
reemplazando en las ecuaciones (37) y (39) v/* 
y ul por pe 

El sistema (35) forma una matriz tridiago- 
nal y puede ser resuelto de manera eficiente 
mediante el algoritmo de Thomas (ver 
Zataráin et al., 1998), una vez especificadas las 
condiciones inicial y de frontera. Es necesario 
señalar que k, en la discretización de la 
condición de radiación fractal, depende de la 
propia solución; sin embargo como el proceso 
de solución del sistema (35) es iterativo, este 
parámetro se calcula en función del estima- 
dor precedente. 

De acuerdo con Zataráin et al. (1998), la 
discretización del dominio se realiza de modo 
= Óx sea constante 
para ¡ =4, 5... N- 2, excepto en la vecindad de 
los drenes; es decir, para x, =0: a) x,-— x, =0.4 
Ox, x,—x,=0.6 Sx, Ax, = 0.1 9x, Ax, =0.6 Ox; y b) 
Xy = 0.4 Ox, Ax,,, = 0.6 Óx, Ax, = 0.1 
9x. El valor de interpolación en el espacio se 


que el incremento x, — X,, 


Xy =L, Xx, 
toma como y= Y en el dominio, excepto, como 
se puede inferir, en la primera y última celdas. 

En cuanto a la discretización del tiempo, 
dada la del espacio, se sigue el enfoque clásico 
de escribir las ecuaciones del movimiento 
en forma adimensional, válido en medios 
homogéneos, para obtener relaciones entre la 
escalas espaciales y temporales características. 
Introduciendo variables adimensionales en la 
ecuación de Boussinesq, ecuación (3), definidas 
como, = XL, <=, H, = A/A, =pfo,R.. 
= R LYTH, donde v, = v(H) y T.= KH, se 
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obtiene la misma ecuación de Boussinesq con 
variables con asteriscos si 1 = V LIT. Debido 
a la naturaleza parabólica de la ecuación 
diferencial, se define el parámetro M = (Ax,)?/ 
At., que puede ser encontrado comparando la 
solución en diferencias finitas con soluciones 
analíticas. El valor del parámetro para los 
tiempos cortos recomendado por Zataráin et al. 
(1998) es del orden de M= 0.1. 


Comparación con una solución analítica 


Con la finalidad de definir valores primeros de 
los parámetros de interpolación en el espacio 
y en el tiempo (y y 0), la solución numérica se 
compara con una solución analítica obtenida 
de la ecuación de Boussinesq en un caso 
particular. Esta solución ha sido construida 
para una linealización de la ecuación diferen- 
cial representada por una transmisibilidad 
constante, pero con una condición de radiación 
lineal en los drenes y que incluye la ecuación 
clásica de Glover-Dumm (Dumm, 1954). Los 
valores utilizados para la simulación son los 
reportados por Fragoza et al. (2003): L =50 m, 
K_=0.557 y m/d, E = 0.1087 m?*/m?, T =2.5065 
m/d,D.=35m,H =50myx=1.,5, 


En la figura 3 se muestra el comportamiento 
de la evolución del abatimiento de la superficie 
libre en un día con diferentes valores de 
interpolación en el tiempo (0). Se muestra 
la variación del abatimiento de la superficie 
libre en todo el dominio de solución y un 
acercamiento sobre el dren, con diferentes 
pasos de interpolación. Puede apreciarse que 
el paso de interpolación óptimo que hace que 
la solución numérica coincida con la solución 
analítica, dado un criterio de error, es (1) = 0.95. 
En la figura 4 se muestra el abatimiento de la 
superficie libre y la evolución del volumen 
drenado por unidad de área de suelo. Los 
resultados muestran que no existen diferencias 
significativas entra la solución analítica y la 
solución en diferencias finitas; se ha utilizado 
M=0.04. 


Comparación de los dos esquemas numéricos 


Los esquemas mixto y en carga se comparan 
entre sí, aceptando los valores y = 0.5 y O) = 
0.95. El suelo utilizado es el caracterizado por 
Saucedo et al. (2003), con los valores 0 = 0.5245 
cm*/cm?, O, = 0 cm*/cm?, K, = 0.446 m/d; los 


valores de los parámetros de las características 
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Figura 3. Evolución de la carga con diferentes pasos de tiempo: a) abatimiento de la superficie libre en un día, 


b) abatimiento de la superficie sobre el dren en un día. 
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Figura 4. Comparación entre la solución numérica y la solución analítica: a) abatimiento de la superficie libre, 


b) evolución de la lámina drenada acumulada. 
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Figura 5. Comparación de los dos esquemas numéricos: a) abatimiento de la superficie libre; 
b) evolución de la lámina drenada, usando las características hidrodinámicas de Fujita-Parlange. 
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Figura 6. Comparación de los dos esquemas numéricos: a) abatimiento de la superficie libre; 
b) evolución de la lámina drenada, usando las características hidrodinámicas de van Genuchten. 
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hidrodinámicas son: a) para Fujita y Parlange 
1, =0.521 m y a. = 0.98; b) para van Genuchten, 
con la restricción de Burdine (1953) m = 1 — 
2/n, m = 0.066 y Y, = -0.15 m. Para comparar 
los esquemas, se propone un distanciamiento 
entre drenes L = 25 m y una profundidad de 
drenes H,=1.5 m. 

Los resultados de la simulación numérica 
obtenidas con las dos características 
hidrodinámicas se muestran en las figuras 
5 y 6. En la figura 5 se tiene la evolución del 
abatimiento de la carga y la lámina drenada 
para tiempos de 60 y 250 d, respectivamente, 
utilizando las características hidrodinámicas 
de Fujita y Parlange, y en la figura 6 se aprecian 
los resultados obtenidos con las características 
hidrodinámicas de van Genuchten para los 
tiempos mencionados con anterioridad. Puede 
verse que no existen diferencias aparentes 
entre el esquema mixto y el esquema en carga 
con cada una de las características hidro- 
dinámicas utilizadas; los resultados muestran 
que utilizar cualquiera de los dos esquemas 
conduce al mismo resultado. 


Problema inverso 
Para evaluar la capacidad de la condición 


de radiación fractal con capacidad de alma- 
cenamiento variable se hace uso de la 


información experimental presentada por 
Zavala et al. (2003). La prueba se realizó en 
un contenedor rectangular con las siguientes 
características: profundidad de los drenes H_= 
120 cm; profundidad del estrato impermeable 
D, = 25 cm, separación entre drenes L = 100 
cm; diámetro y longitud de los drenes d = 5 
cm y | =30 cm. Las propiedades del suelo son 
porosidad volumétrica q = 0.539 cm?*/ cm? y 
conductividad hidráulica saturada K, = 18.3 
cm/h. La dimensión fractal relativa del suelo 
obtenida con la ecuación es s = 0.7026. 

Las características hidrodinámicas utili- 
zadas son las de Fujita y Parlange; para esto, 
se fija el valor de a. = 0.98 y los parámetros 
A, y K, se estiman mediante la minimización 
de la suma de los cuadrados de los errores 
entre la lámina drenada medida y la lámina 
drenada calculada con la solución numérica 
en el transcurso del tiempo. En la figura 7 
se presentan los resultados obtenidos del 
problema inverso para 24 h y 240 h. Se usaron 
los dos esquemas numéricos y el resultado fue 
el mismo. La mejor aproximación a los datos 
experimentales es obtenida con A, = 0.5152 m 
y k,= 0.10, que proporcionan un ECM = 0.1283 
cm. 

La condición de radiación fractal y la capa- 
cidad de almacenamiento variable reproducen 
de buena manera los datos medidos en labo- 
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Figura 7. Evoluciones de la lámina experimental y calculada con la condición de radiación fractal 


y capacidad de almacenamiento variable, utilizando el esquema numérico propuesto. 
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ratorio para el intervalo de tiempo menor a 6h, 
posteriormente hay una ligera subestimación 
de los mismos hasta un tiempo de 150 h; sin 
embargo, el buen acuerdo entre la lámina 
drenada medida y la lámina drenada obtenida 
con el modelo es evidente. 


Conclusiones 


Se ha resuelto la ecuación unidimensional de 
Boussinesq del drenaje agrícola con el método 
de diferencias finitas basada en un balance 
local de masa. Resultaron dos esquemas de 
discretización de la derivada temporal, la cual 
representa el cambio de almacenamiento en este 
balance. En uno aparecen de manera explícita 
la carga y porosidad drenable, variables 
ligadas con una relación funcional, que se ha 
denominado esquema mixto; en el otro aparece 
sólo la carga hidráulica, denominado esquema 
en carga. Los dos esquemas coinciden cuando 
la porosidad drenable es independiente de la 
carga. 

La validación parcial de ambos esquemas 
fue realizada mediante la comparación de 
las soluciones numéricas obtenidas con una 
solución analítica de la literatura construida 
para condiciones de linealidad. Las evolu- 
ciones de la carga de agua en el perfil y la 
lámina drenada calculadas con la solución 
analítica son similares, bajo un criterio 
de error, a las calculadas con la solución 
numérica para todo tiempo. La aplicabilidad 
de los esquemas propuestos para condiciones 
donde la linealidad es restringida, la ausencia 
de fluctuaciones tanto en el tiempo como en 
el espacio de la carga y de la lámina drenada 
permite recomendar los esquemas numéricos 
de la ecuación unidimensional de Boussinesq 
propuestos para el estudio de la dinámica 
del agua en los sistemas de drenaje agrícola 
subterráneos. En particular la solución 
numérica construida puede ser utilizada para 
la caracterización hidrodinámica del suelo a 
través de una modelación inversa, es decir, 
a partir de las evoluciones experimentales se 
pueden inferir los parámetros del sistema. 


La solución numérica propuesta puede ser 
utilizada para un mejor diseño de sistemas 
de drenaje agrícola subterráneo, ya que las 


hipótesis consideradas en las soluciones 


clásicas han sido eliminadas. 
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Abstract 


CHÁVEZ, C., FUENTES, C. € ZAVALA, M. Finite difference solution of the agricultural 
drainage Boussinesq equation with variable drainable porosity subject to a fractal radiation 
boundary condition. Water Technology and Sciences, formerly Hydraulic engineering in 
Mexico (in Spanish). Vol. I, No. 4, October-December, 2010, pp. 105-117. 


Subsurface drainage systems are used to control the depth of the water table and to reduce or 
prevent soil salinity. Generally, the flow of the groundwater is studied with the Bousssinesq 
equation, whose analytical solutions are obtained assuming that aquifer transmissivity and 
drainable porosity are constant. These solutions assume as well that the free surface of the water 
falls instantly over the drains. The general solution requires numerical methods. Some authors 
have demonstrated that the drain boundary condition is a fractal radiation condition and that 
the drainable porosity is a variable which is related to the soil retention curve. This solution 
has been obtained with a finite element method, which in one-dimensional form is equivalent 
to a finite difference method. Here, we propose a finite difference solution of the differential 
equation with variable drainable porosity and a fractal radiation condition. The proposed finite 
differences method has two formulations: the first one, with an explicit head and drainable 
porosity, both joined with a functional relationship, which we call mixed formulation; and 
the second one, which we call head formulation, with only the head. Both methods have been 
validated with a lineal analytical solution, and the nonlinear part is stable and brief. The 
proposed numerical solution is useful for the hydraulic characterization of soils with inverse 
modeling and for improving the design of agricultural drainage systems, considering that the 
assumptions of the classical solution have been eliminated. 


Keywords: mixed formulation, head formulation, soil water retention curve, inverse 
modeling. 
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